Synchronization in electric power networks with inherent heterogeneity up to 100% inverter-based renewable generation

The synchronized operation of power generators is the foundation of electric power network stability and a key to the prevention of undesired power outages and blackouts. Here, we derive the conditions that guarantee synchronization in power networks with inherent generator heterogeneity when subjected to small perturbations, and perform a parametric sensitivity analysis to understand synchronization with varied types of generators. As inverter-based resources, which are the primary interfacing technology for many renewable sources of energy, have supplanted synchronous generators in ever growing numbers, the center of attention on associated integration challenges have resided primarily on the role of declining system inertia. Our results instead highlight the critical role of generator damping in achieving a stable state of synchronization. Additionally, we report the feasibility of operating interconnected electric grids with up to 100% power contribution from inverter-based renewable generation technologies. Our study has important implications as it sets the basis for the development of advanced control architectures and grid optimization methods that ensure synchronization and further pave the path towards the decarbonization of the electric power sector.

The simultaneous differential equations that describe the motion of the masses in the system depicted in Supplementary Figure 1, assuming initial values equivalent to zero, are: Assuming F = 0, in order to study the natural modes, they can be rewritten as In matrix form, defining x = x 1 x 2 T , they are: By defining v 1 =ẋ 1 and v 2 =ẋ 2 , a first order state space model is created:ẋ Solving this system of differential equations yields the dynamical equations governing the speed as a function of time, v 1 (t) and v 2 (t), respectively. Then kinetic energy as a function of time can be computed by E(t) = 1 2 m 1 · v 2 1 (t) + 1 2 m 2 · v 2 2 (t). Supplementary Table I

Supplementary Note 2: Dynamic model of the generators
Electric power systems can be modeled as a complex network with power generators, substations, and load centers constituting the nodes and the transmission lines that interconnect them the links. Generally, the behavior of power systems can be described using a set of differential-algebraic equations (DAE) equations, as in [1]:ẋ (t) = f (t, x(t), y(t), u(t)) 0 = g(t, x(t), y(t), u(t)) (5) where dot notation indicates time derivative;ẋ = d dt x, t denotes time, x is the vector of state variables, y is the vector of algebraic variables, and u is the vector of inputs. The vector field f (t, x, y, u) describe the dynamic evolution of the state variables of the system including generation, load, and any other dynamic elements and the associated control systems. The algebraic equations determined by g(t, x, y, u) describe the system constraints and the electric network that interconnects the dynamic elements [1]. In this study, we focus on the impacts of generator parameters on the system synchronization. Thus, we model the system generators [1] as dynamic elements. The transmission lines and load centers, whose contribution to inertia is relatively minimal, are modelled as algebraic elements. The power generation technologies considered here can be primarily divided into two main categories: (1) synchronous generator-based power plants and (2) power electronics-based power plants.
Synchronous Machine-based Power Plants: In conventional power systems, the majority of generation facilities are equipped with synchronous generators that transform mechanical power directly into alternating current electric power. Over the past century, synchronous generators -including various types driven by disparate fuels such as coal, natural gas, nuclear fission, hydro, etc -have matured and presently offer a reliable and economical performance for the bulk generation of electricity. They have two main electromagnetic parts: (1) Stator: which is a set of stationary windings, sinusoidal distributed around a rotating shaft and (2) Rotor: which is a mechanically rotating shaft housed within the stator, that has the primary field winding, and damper windings. The damper windings on the rotor shaft help to damp out the machines oscillations induced during transient behaviors [2]. The general structure of a synchronous generator is shown in Supplementary Figure 2. When the loading condition on synchronous generator changes, the rotation speed, and thus the frequency changes. This change in frequency is registered by a governor, which varies mechanical power to match the electric power demand, through a feedback control process with the objective to track the frequency reference. The adjustment of mechanical power comes in three steps. First, the mechanical power adjustment is determined by the droop value, which is a linear compensator. Second, the mechanical power is adjusted accordingly by the governor, which elicits a turbine response. Third, the frequency excursion is arrested when the adjusted mechanical power applies to the shaft of the synchronous machine. This process, assuming the voltage dynamics remain ideal, is summarized in Supplementary Figure 3.
Using simplification rules for closed-loop feedback control systems, one can write a transfer function of the synchronous generator frequency response to changes of power output as: This transfer function presents a reduced-order model for the frequency response of the synchronous generator. It should be noted that D m is a function of rotor angle (δ), D m = f (δ), though for small deviations from the synchronous speed, its variation is minimal and can be assumed a constant, D m (δ) ≈ D m [2], [3].
The roots of the characteristic equation for the synchronous generator are calculated by Looking at (s 1,2 ), the damping component is mainly determined by the real part of these roots, −( 1 T T G + Dm R·Mm ). According to this relationship, for a synchronous generator, the components that contribute to damping are (1) the mechanical damping provided by the damper winding, (2) product of droop gain and mechanical inertial momentum, and (3) the response time constant of the governor and turbine system. Also from the transfer function, one can draw a conclusion that the pace of frequency response (the inertial component) is directly tied to the mechanical inertia ( 1 Mm ), particularly the fast transient portion of the response because of the existence of a zero (the term s in the numerator). The mathematical description of the frequency response of a synchronous generator concludes the followings: the larger the mechanical dampers, D m , and/or the smaller the product of the droop gain and mechanical inertia, (R.M m ), and/or the smaller governor and turbine response time, T T G , (indicating a faster response time), the larger the magnitude of the negative real part of imaginary number, which indicates a higher damping of the natural modes. The larger the mechanical inertia, (M m ), the slower the rate of change of frequency, a well known result in power system dynamics.
We verified and validated this reduced-order model for synchronous generator frequency response, described in Eq. (6), using higher-order models developed in PSCAD (abbreviation of "power system computer aided design"). We used the open-source models developed by the National Renewable Energy Laboratory (NREL) and University of Colorado Boulder (CU Boulder), which are available to public for research use [4]. The mathematical description of the high-order models for synchronous generators and associated power network used in these PSCAD test cases are explained in [5]. PSCAD is an industry grade software package that operates in the electromagnetic transient (EMT) domain and is capable of capturing the fast dynamics of controllers that are generally obscured in phasor domain simulation tools. All three electrical phases are modelled with a suitable time resolution (in the microsecond range) to capture fast transients. The results of the validation are shown in Supplementary Figure 4 and corroborate that the reduced-order model is effective at approximating the frequency response of a synchronous generator. 59 Figure 4: Validation of the synchronous generator frequency response using the reduced-order model by comparing it with the frequency response using the high-order results obtained from the PSCAD simulation for 3.5MW load switching in the 9-bus system -the model is available open-source at [5].
The transfer function of the reduced-order frequency response model is a second-order dynamical system. Therefore, for studying synchronization, we use a generic second-order model in the form of a second-order ordinary differential equation (common harmonic oscillator model) with electric angle and speed as the dynamic variables. Using the generic model, assuming x 1 = δ sg , x 2 = ω sg , the frequency response of a synchronous generator can be described by: where p * and p e are the mechanical power input and electrical power export, respectively. The M sg and D sg are inertia and damping coefficients.These coefficients represent all contributing elements to inertia and damping components. Because of the mechanical nature of energy conversion with synchronous machines, the response agility is bounded by relatively slow mechanical adjustments. The inertia and damping coefficients are constant, as they are dependent on the machine design, and can be expressed by the D/M ratio for simplicity. The values of these variables are physical constraints and are set by the machine and power plant design; during system operation, they generally cannot be changed. When subjected to a disturbance that results in mismatch between the generator mechanical power and electric power output, the electric speed and subsequently the electric angle, experience transients. We verified and validated this generic second-order model presented in Eq. 8 using the previously validated reduced-order model presented in Eq. (6). The results of the validation are shown in Supplementary Figure 5 and corroborate that the generic second-order model can effectively approximate the frequency response of the synchronous generator. 59

Reduced-order Model Generic Second-order Model
Supplementary Figure 5: Validation of the synchronous generator frequency response obtained from the generic second-order model by comparing it with the frequency response obtained from the reduced-order model, emulating a 3.5MW load switching event in the 9-bus system .
Power electronics-based power plants: Power electronic inverters, which are often used to interface variable renewable energy sources with the grid, can be categorized by the primary control mode as either (i) grid-following (GFL) or (ii) grid-forming (GFM). Both classes of inverters have similar hardware components with the structure differing chiefly in the control schemes, which has a considerable impact on the interactions with the electric power grid. GFL inverters, the most commonly used class of inverters, rely on the phase locked loop (PLL) to estimate and track the frequency at the point of interconnection (POI). Typically, the control objective of these devices is to maintain a set power export regardless of varying grid conditions. If operated as a grid support device using a droop-relationship, then active power and reactive power are adjusted proportionally to frequency and voltage deviations. In contrast, GFM inverters, an emerging and promising technology as applied to the parallel operation in bulk power systems, adjust the frequency and voltage according to the power injection deviations at the POI. In other words, they are capable of supporting the grid by generating frequency independent of the grid operating conditions. An additional distinction in GFM control is the necessity of available positive headroom to perform these control objectives. Although there are variants for each of these two classes in the literature [6]- [9], this work focuses on the two most common state-of-the-art classes: (1) the GFM with either the virtual-inertia support functionality or the multi-loop droop control approach and (2) the GFL. The general structure of these technologies is shown in Supplementary Figure 6.
The GFM with virtual-inertia support functionality is also known as the virtual synchronous machine (VSM) [10], with the aim to emulate a synchronous generator frequency response to mitigate Figure 6: Comparative illustration of simplified GFL and GFM control systems potential operational challenges that grid-coupled inverters may pose [11]- [18]. Subsequently, at least 11 control approaches have been proposed for the VSM [10], [11], [19]- [38], all of which are neatly summarized in [39], with the primary differences generally associated with internal decision-making and the choice of state variables. Nonetheless, all of these various approaches have a similar output as they try to emulate a synchronous generator and, thus, their frequency response can be described by the second-order differential equation used for a synchronous generator. The main distinction is that the physical variables should be replaced with their digital dual and a similar frequency response mechanism applies.
The mechanism of frequency response in the multi-droop GFM inverter adjusts the frequency set point as a function of output power deviations through an open-loop control process that leverages the frequency-droop relationship of real power flow across AC transmission systems. The adjustment of frequency comes as the following: Output power is determined by the voltage values at the two terminating points of the grid coupling inductor of the LCL filter. Denoting E and V the two voltage values and the value of the inductance X Lc , the power transfer is P g = E·V X Lc sinδ. Assuming the voltages are very close to 1 p.u. and after linearization, the change of the multi-loop GFM power output can be described by ∆P g = 1 X Lc ∆δcosδ. Since δ is a small value, cosδ ≈ 1, and ∆P g = 1 X Lc ∆δ. For frequency response, first, the power error is detected through a low-pass filter whose output is used to create the power error signal for the droop controller. Second, the droop controller determines the output frequency set point. Finally, passing the frequency set point through an integrator creates the angle of voltage at the capacitor of the LCL filter. As noted, this is an open-loop process. This process is summarized in Supplementary Figure 7.
One can write a transfer function of the multi-loop GFM frequency response as a function to changes of power demand as: where • R is the droop gain, and • T ωc is the low-pass filter time constant, the reciprocal value of cut-off frequency (ω c ), The characteristic equation for the multi-loop GFM inverter is a first-order polynomial whose singular root is which is the reciprocal value of the time constant of the low-pass power measurement filter, equal to the cut-off frequency. Knowing the cut-off frequency is relatively high Tω c > 100 rad/sec) [40], the only root of this equation is well to the left hand side of the imaginary axis and always real, and thus, the multi-loop GFM inverter is safe be assumed asymptotically stable.
The transfer function indicates that the frequency response of multi-loop GFM is first-order. The pace of response (the inertial component) is directly determined by the cut-off frequency of the lowpass power measurement filter, ω c , and the droop gain R. The settling value is determined by the droop gain, R; the higher the cut-off frequency or larger the droop gain, the faster the frequency deviates, and the smaller the droop gain, the smaller the frequency deviation.
We verified and validated this reduced-order model for GFM frequency response, described in Eq. (9), using the high-order models developed in PSCAD. We used the very same open-source models developed by the NREL and University of Colorado Boulder (CU Boulder) [4] that previously were used to verify the synchronous generator reduced-order model. The mathematical description of the high-order models for the GFM inverter and the associated power network used in this PSCAD test case are detailed in [5]. The results of validation are shown in Supplementary Figure 8 and indicate that the reduced-order model can suitably approximate the frequency response of the GFM.
The transfer function of the reduced-order frequency response model for the multi-loop droop GFM yields a first-order system which can be approximated by order-reduction of a second-order system, while the transfer function for the VSM-GFM yields a second order model. Therefore, for studying synchronization, we use a generic second-order model in the form of a second-order ordinary differential equation (common harmonic oscillator model) with electric angle and speed as the dynamic variables. Using the generic model, assuming x 1 = δ pe , x 2 = ω pe , the frequency response of GFM can 59 Figure 8: Validation of the GFM frequency response obtained from using the reducedorder model obtained from the analytical model by comparing it with the frequency response using the full-order results obtained from the PSCAD EMT simulation for a 3.5MW load switching event in the 9-bus system be described byẋ where p * and p e are the power output setpoint and actual power export, respectively. The M pe and D pe are the equivalent inertia and damping coefficients; when {M >> 0} the VSM control is employed and for {M ≈ 0 ∧ M ̸ = 0} the control mode is the multi-loop droop class of GFM (which would result in order reduction to the first-order). These coefficients represent all contributing elements to inertial and damping component of frequency response In power electronics-backed generation units, because of the non-mechanical nature of energy conversion with DC sources, the response agility is primarily bounded by the digital controller process time, which could be very quick, and can be modified in an adaptive manner and in real-time, except for the physical value of the grid-coupling inductance. Therefore, unlike synchronous generators, the inertia and damping coefficients in GFM are independent of each other. The power setpoint defines the amount of power that is commanded to be released through the semi-conductor switches. This model implies that following a disturbance that results in mismatch between the GFM power setpoint and power output, electric angle can begin to oscillate. In a stable operation, these oscillations should be damped out within a relatively short time following their inception -within seconds. We verified and validated this generic second-order model presented in Eq. 11 for multi-loop GFM using the previously validated reduced-order model presented in Eq. (6). The synchronous generator model validation accounts for the VSM-GFM model validation. The results of the validation are shown in Supplementary Figure 9 and show that the generic second-order model can effectively approximate the frequency response of a multi-loop GFM.
The GFL, on the other hand, does not contribute to the system electric angle construction. The GFL does not produce any constituting voltage frequency and instead follows the frequency at the POI. It estimates the frequency at the POI using a PLL and through a linear relationship (droop) it adjusts its power export by acting as a current source. Therefore, the its governing equation can 59

Reduced-order Model Generic Second-order Model
Supplementary Figure 9: Validation of the multi-loop GFM frequency response obtained from the generic second-order model by comparing it with the frequency response obtained from the reducedorder model, emulating 3.5MW load switching in the 9-bus system be given by: where p * and p e are the power setpoint and power export. Effectively, the GFL can be modeled as a negative load.

Supplementary Note 3: Linear model of electric power network
To obtain the Jacobian of the electric power network of interest in this study, first, let us consider the operation at a stable equilibrium where the generator power setpoints and electric power exports are equal, P * = P e , for both the GFM and Synchronous generators. Following a disturbance, there will be a mismatch between these known as the acceleration power, defined by P = P * − P e as pertaining to changes in generators' power output. Inspired by Machowski's [2] analysis of power networks with only synchronous generators considered, linearization of the governing equations for the GFM and Synchronous machines (the linearization introduces the linear difference operator, ∆) interconnected at the i-th node (8) and (11) yields: where ∆P i is the change in the generator power injection. For any node i, the active power injection at the node with voltage of V i = |V i |∠δ i can be expressed by power-flow equation as [2]: where Y ij = G ij + j B ij here are the effective admittance values obtained from Kron reduction; G is the conductance and B is the susceptance. The linearized form of (14) can be expressed in a matrix Equation (15) can alternatively be presented as: 1 H 1 ∆P n = H n ∆δ n (17) where H is the system network Laplacian (effective admittance), an (n × n) matrix whose elements can be computed by: Replacing (17) into (13)  In (19), I is the n × n identity matrix and 0 represents a (n × n) matrix of zeros. This is a (2n) set of equations representing n electric angles and electric speeds of all n generators.

Supplementary Note 4: Coordinate transformation
In power system stability analysis, one bus is appointed as the reference generator (commonly, the slack generator) and all other angles are analyzed relative to the angle of the reference generator [2], [3], [41]. Subsequently, system stability is analyzed as relative angle stability that assesses whether the angles remain relatively bounded, which does not hold true if the pace of their evolution is not homogeneous (i.e., the definition of synchronization). This is a common technique and here we borrow our analysis from [2]. From (18) it can be seen that H is singular: This singularity can be used for a coordinate transformation to express all electric angles as relative electric angles with respect to the reference angle. Let us suppose that the reference angle is the electric angle of the n-th generator. For the coordinate transformation, one can write [2]: In (21), for the last term, H in ∆δ n , by leveraging the singularity property of the Laplacian it can be expressed as H in ∆δ n = − n j=1,j̸ =n H ij ∆δ n . Substituting it in (21) produces Now, (15) can be expressed in the form of relative angles as: And in compact form: where H n−1 is an n × (n − 1) matrix and the angles presented in ∆δ i,n are the relative angles between node i and the reference angle at the reference generator n. Replacing (24) into (13) produces: In (25), [0] represents a (n × n) matrix of zeros and [1 −1 ] represents a matrix of identity with an extended column whose arrays are all equal to (-1). This is a (2n − 1) set of equations representing (n − 1) electric angle of all generators relative to the reference generator and n actual electric speed that would include the reference generator's electric speed.

Supplementary Note 5: Deriving the stability condition
Here we solve the necessary condition for stable synchronization in electric power networks with heterogeneous coupling damping factors. The linear model previously developed and described in (25) can alternatively be described as: where ∆ is the linear difference operator, δ i,n is the vector of (n − 1) relative electric angles and ω i and ω n are absolute electric speeds expressed by a vector of (n − 1) and an integer value (making it a vector with a single array), respectively. [0] is a matrix of zeros, and [I] is the identity matrix, each a square matrix of dimension (n − 1), and [−1] is a matrix with dimension 1 × (n − 1) and all elements are (-1).
each a square matrix with dimension of (n − 1), with H, D, and M the network interconnection Laplacian coefficients and generators damping and inertia coefficients, respectively. Finally, h n = −M −1 n H i,n , ∀i = 1, 2, · · · , (n − 1) is a 1 × (n − 1) matrix and d n = −M −1 n D n is a single array.
By definition of eigenvalues: Jw = λw (28) with J being the system Jacobian along the solution trajectory. It can be extended as This can be presented in compact form as Expansion of this expression provides: where · is the inner product operator and 1 is a matrix of all ones of (n − 1) × 1.
From the third equation in (31), it can be derived: From the first equation in (31): Replacing (32) into (33) produces where ⊗ is the outer product operator. Let us rearrange the second equation from (31): Replacing (34) in this expression produces: Since w ′ cannot be zero, then the first expression should be equal to zero. Expression (36) can be expanded to: where λ are the eigenvalues of the system. The necessary condition for existence of eigenvalues is that this characteristic matrix is singular and thus, det(p(λ)) = 0. The existence of stable equilibria is guaranteed if the all the eigenvalues have non-positive real part. One should note that although p(λ) is a third-order polynomial, the solution yields only a total of (2n − 1) eigenvalues whose eigenvectors are linearly independent for a network with n generators. This is because Eq. (32) of the above procedure projects the (2n − 1)th eigenvector into 2(n − 1) other dimensions, creating a total of 3(n − 1) eigenvectors. But at the end, there are only (2n − 1) eigenvalues whose eigenvectors are linearly independent and the others are repetitive of the (2n − 1)th eigenvalue and the corresponding eigenvector.

Supplementary Note 6: Analytical derivation of stability condition for 3-generator system
Let us consider a 3-generator system (classic 9-bus system), depicted in Supplementary Figure  10. This is a modified system to adopt the reduced-order generator models, developed based on a standardized benchmark without any change to the line parameters or network topology [3]. Supplementary Figure 10: One-line diagram of the 3-generator, 9-bus system, and the dynamic and network data For this system, n = 3 and, let's assume gen3 as the reference generator. By definition of eigenvalues, using its Jacobian, it can be said: From the last equation in (39): From the first two equations in (39), it can be derived: Replacing (40) into (41) yields: From the third and fourth equation in (39): Replacing (42) in (43) produces: This expression can be rearranged as: It can be expanded as: And rearranging produces: In compact form, it can be written as The eigenvalues of system can be found by solving the det(ψ) = 0 for λ. We have demonstrated the results for this test system using three different combinations of generation technologies and have validated them using the industry-grade full-order model, shown in Table I of the Main Manuscript. The coefficients for the reduced-order frequency response models corresponding to these case and subsequent nadir frequency and rate of change of frequency (ROCOF) are presented in Supplementary Table II Recalling Eq. (13), let us suppose a special condition where the damping factors, d i = − D i M i , are homogeneous, constituting d 1 = d 2 = · · · = d n . The stability criteria for this special condition is developed by Motter [42] and Machowski [2]. For the sake of argument, we borrow their analysis here. The assumption of homogeneity leads to a reduced system Jacobian to 2(n − 1) by representing electric speeds as relative electric speeds with respect to the electric speed of the reference generator. To this end, let us assume the coupling damping factor is uniform for all generators across the system and equivalent to (− D M ). For any i-th generator with the n-th generator being the reference generator, it can be written [2]: where ω i,n = ω i − ω n is the electric speed deviation from the reference machine. Substituting (49) into (13) produces where ∆δ n−1 and ∆ω n−1 are vectors of (n − 1) relative electric angle and electric speeds and [0] and [I] are (n − 1) × (n − 1) matrices of zero and identity. d is a diagonal matrix of damping coefficients with all of its diagonal arrays equivalent to (− D M ) and h is a (n − 1) × (n − 1) matrix of relative network coefficients.
To assess the stability of this system, one can calculate the eigenvalues of (50), denoted by λ with corresponding eigenvectors of w. The definition of eigenvalues and eigenvectors provides [2], [42]: is the eigenvalue of the submatrix h. It can be rearranged as Equation (50) is equivalent to (26) under the assumption of d i homogeneity yielding a zero coupling mode k d i = 0 in (50). The eigenvalues of (50) can be found by solving (52) for λ as

Supplementary Note 8: Frequency response in power networks and hertz-sec metric
The power network operator is concerned with continuously maintaining a stable and homogeneous system frequency across the network, 60Hz in North America and 50Hz in the rest of the world, reflecting the balance between the electric power generation and consumption. In real-world instances, this system frequency is continuously fluctuating within a small, stable equilibria as a result of slight, countless perturbations, but can be approximated as a quasi-steady state system frequency due to the small boundedness of the deviations. Any perturbation can induce transient oscillations and in a stable system these oscillations are damped out within a relatively short time following their inception. As an example, the frequency response trace from the European interconnected grid is shown in Supplementary Figure 11. The largest value of frequency excursion is called the nadir, which is an important value in power grid management and control because if the frequency drops below a standard threshold -48.5Hz in Europe [43] and 59.5Hz in North America [44] -it can trigger protective equipment that results in an interruption of power supply to some parts of the grid as a preventative measure to avoid cascading failures and blackout. Another important frequency related metric is the rate of change of frequency (ROCOF), where high values can trigger protective equipment. This metric quantifies the tangent of frequency response right after disturbance occurs and, subsequently, when the frequency starts to change. The ROCOF is measured in Hz sec . The frequency nadir and ROCOF play a critical role in the operation of power networks and, but the individual values do not provide sufficient information about the frequency response and dynamic state of system. Therefore, to quantify and assess the frequency response, we use the hertz-sec metric [45] which is a proxy for the changes of kinetic energy induced by the disturbance. This metric integrates the absolute value of frequency deviation over its transient period. The hertz-sec, HS, is defined as: where f 0 and f (t) are the pre-disturbance frequency (a constant) and the post-disturbance frequency (a function of time), respectively, and t 0 and t s are the time the frequency departs the pre-disturbance steady-state value and the time it permanently reaches the settling steady-state value, respectively. A smaller value of this hertz-sec metric, as a proxy for reflecting smaller amount of kinetic energy, indicates a faster and more agile response, which is desired.  Figure 11: Frequency response dynamic for the European synchronized continental power grid recorded by the Phasor Measurement Units (PMUs) in Italy, Belgium, Germany, and France on Jan 8, 2021: In power grid operation, the nadir frequency is the lowest point of the frequency dynamics and it is a critical value because protective equipment are set to operate if the nadir value drops below a certain threshold to avoid cascading failures and a potential major blackout. To prevent loading-shedding or disconnection of part of the grid, it is important that at dynamic state of the grid is maintained such that all times, the nadir frequency for most common small disturbances remains above this threshold level.

Supplementary Note 9: Power network data
We used the network data from the standard IEEE benchmarks including IEEE 9-bus [3], IEEE 30-bus [46], IEEE 39-bus [47], IEEE 57-bus [48], IEEE 118-bus [49], and IEEE 300-bus [50] systems. The steady state model for these systems is freely available as part of the MATPOWER tool [51] and a statistical description of these systems is provided in Supplementary Table III. Supplementary The dynamic data for some of these benchmarks does not exist. For this study, we generated synthetic dynamic data which represent real-world power system dynamics and characteristics relying on the data available from those publicly available benchmarks [3], [52]- [54]. Creation of synthetic power grids is a common practice because the real-world power system data is not always easily accessible due to legal matters, pursuant to Homeland Security Act of 2002: Critical Infrastructure Information Act [55]. Variety of dynamic power systems have been developed and documented [56]- [61], mostly extremely large systems which have application for optimization problems and power flow analysis. Our synthetic dynamic data are generated algorithmically in base cases to represent conventional power network with all generation stations equipped with synchronous generator. The range we chose is consistent with those empirically observed in the data available. The total power (MVA) capacity of generators are also randomly allocated between 100% and 250% of their total power dispatch determined by the power-flow equations to account for reserve margin. Then we validated by perturbing all generators, one generator at a time, and monitored the dynamic response of all generators that are interconnected to the network to ensure they are stable and their characteristics and response match what is commonly observed in power systems. To perturb the system and emulate the load switching, we subjected the electric angle of generators, δ to a step change that replicates a change of the loading condition and subsequently, we monitored the electric speed, ω. Because we applied Kron reduction which is a common technique in power network synchronization studies, the load nodes are eliminated and there is no other mean to directly switch loads, other than switching the δ of generators in a sustained fashion. System frequency can be measured as the electric speed of generators, ∆ω = ∆δ, as ∆f = ∆ω(2π) −1 + f ss with f ss = 60Hz (in the United States, 50 Hz in Europe), and (∆δ 0 , ∆ω 0 ) = (0, 0), assuming the frequency departs from an equilibrium when disturbed.